
import nbconvert


jupyter nbconvert --to html HW2_.ipynb

  File "<ipython-input-3-c97678e9646a>", line 1
    jupyter nbconvert --to html HW2_.ipynb
                    ^
SyntaxError: invalid syntax


from OSMPythonTools.nominatim import Nominatim
import overpy
import shapely.geometry as geometry
from shapely.ops import linemerge
import folium
import random


nominatim = Nominatim()
areaId = nominatim.query('Afganistan').areaId()
areaId

3600303427


query = """[out:json][timeout:25];
rel(303427);
out body;
>;
out skel qt; """

api = overpy.Overpass()
result = api.query(query)

lss = [] #convert ways to linstrings

for ii_w, way in enumerate(result.ways):
    ls_coords = []

    for node in way.nodes:
        ls_coords.append((node.lon, node.lat)) # create a list of node coordinates

    lss.append(geometry.LineString(ls_coords)) # create a LineString from coords


merged = linemerge([*lss]) # merge LineStrings
afganistan = geometry.Polygon(merged)


way = []
lat = random.uniform(-90, 90)
long = random.uniform(-180, 180)
way.append(geometry.Point(long, lat))


while True:
    if afganistan.contains(geometry.Point(long, lat)):
        break
    if (lat <= afganistan.centroid.y) and (long <= afganistan.centroid.x):
        lat = random.uniform(lat, afganistan.centroid.y)
        long = random.uniform(long, afganistan.centroid.x)
    elif (lat >= afganistan.centroid.y) and (long >= afganistan.centroid.x):
        lat = random.uniform(afganistan.centroid.y, lat)
        long = random.uniform(afganistan.centroid.x, long)
    elif (lat <= afganistan.centroid.y) and (long >= afganistan.centroid.x):
        lat = random.uniform(lat, afganistan.centroid.y)
        long = random.uniform(afganistan.centroid.x, long)
    elif (lat >= afganistan.centroid.y) and (long <= afganistan.centroid.x):
        lat = random.uniform(afganistan.centroid.y, lat)
        long = random.uniform(long, afganistan.centroid.x)
    way.append(geometry.Point(long, lat))


map_of_afganistan = folium.Map([33.83133703361575, 66.02693526354392], zoom_start=2)

for point in range(len(way)):
    folium.Marker(
      location = [way[point].y, way[point].x], 
      popup = str(point) + ' step' 
     ).add_to(map_of_afganistan) 
folium.GeoJson(afganistan).add_to(map_of_afganistan)
map_of_afganistan


from OSMPythonTools.overpass import Overpass
overpass = Overpass()
result = overpass.query('relation["admin_level"="2"][boundary=administrative];out;')
areas = [] 
for el in result.elements():
    el_id = str(el.id()) 
    if el_id[-1] == '0':
        areas.append([el.id(), el.tag('name:en')])


from OSMPythonTools.overpass import overpassQueryBuilder, Overpass
counts = []
for area in areas:
    areaId = 3600000000 + area[0]
    overpass = Overpass()
    query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"landuse"="residential"', out='count')
    result = overpass.query(query, timeout=900)
    counts.append([result.countElements(), area[0], area[1]])


def maximum_residential(counts):
    maximum = counts[0]
    for country in counts:
        if maximum[0] < country[0]:
            maximum = country
    return maximum

max_residential = maximum_residential(counts)
print(f'В {max_residential[2]} (id: {max_residential[1]}) содержится наибольшее количество объектов жилой недвижимости - {max_residential[0]}')
cameroon_id = max_residential[1]

В Cameroon (id: 192830) содержится наибольшее количество объектов жилой недвижимости - 139054


areaId = 3600000000 + cameroon_id
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"amenity"="atm"')
result = overpass.query(query, timeout=900)


points = []
for atm in result.elements():
    points.append((atm.lon(), atm.lat()))


import rtree.index

idx = rtree.index.Rtree()
for i, p in enumerate(points):
    idx.insert(i, p + p, p)

p = (way[-1].x, way[-1].y)
five_nearest = list(idx.nearest(p, 5, objects='raw'))


from shapely.geometry import Point, LineString
orig_point = Point(way[-1])
lines = [LineString([orig_point, Point(dest)]) for dest in points if dest in five_nearest]


from pyproj import Geod
geod = Geod(ellps="WGS84")
lines_length = [geod.geometry_length(line) for line in lines]
print(f'{lines_length} расстояния от точки до 5 ближайших банкоматов в Камеруне')

[6144389.378466479, 6557451.369484478, 6557127.224546507, 6514995.396643748, 6489205.27622738] расстояния от точки до 5 ближайших банкоматов в Камеруне


sum_length = 0
for i in lines_length:
    sum_length += i
print(f'{sum_length} - сумма расстояний в метрах от точки до 5 ближайших банкоматов в Камеруне')

32263168.64536859 - сумма расстояний в метрах от точки до 5 ближайших банкоматов в Камеруне


countries_and_capitals = []
for area in areas:
    areaId = 3600000000 + area[0]
    overpass = Overpass()
    query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"capital"="yes"')
    result = overpass.query(query, timeout=900)
    if result.toJSON()['elements'] == []:
        continue
    else:
        countries_and_capitals.append([area[0],
              area[1],
              result.toJSON()['elements'][0]['id'],
              result.toJSON()['elements'][0]['tags']['name:en']])


countries_and_relation_capitals = []
for country in countries_and_capitals:
    if country[1] == 'Brazil':
        areaId = 3600000000 + 2758138
    elif country[1] == 'Australia':
        areaId = 3600000000 + 2354197
    else:
        nominatim = Nominatim()
        capital_country = country[3] + ', ' + country[1]
        areaId = nominatim.query(capital_country).areaId()
    countries_and_relation_capitals.append([country[0], country[1], areaId, country[3]])


num_pharmacies = []
for country in countries_and_relation_capitals:
    capitalId = country[2]
    if capitalId is None:
        continue
    capitalId = capitalId
    overpass = Overpass()
    query = overpassQueryBuilder(area=capitalId, elementType=['node','way', 'relation'], selector='"amenity"="pharmacy"', out='count')
    result = overpass.query(query, timeout=900)
    num_pharmacies.append([result.countElements(), country])


def maximum_pharmacy(num_pharmacies):
    maximum = num_pharmacies[0]
    for country in num_pharmacies:
        if maximum[0] < country[0]:
            maximum = country
    return maximum

max_pharmacy = maximum_pharmacy(num_pharmacies)
print(f'В {max_pharmacy[1][3]} столице {max_pharmacy[1][1]} содержится наибольшее количество аптек - {max_pharmacy[0]}')
bangladesh_id = max_pharmacy[1][0]

В Dhaka столице Bangladesh содержится наибольшее количество аптек - 1286


areaId = 3600000000 + bangladesh_id
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"amenity"="school"')
result = overpass.query(query, timeout=900)


points = []
for school in result.elements():
    if school.lon() is None or school.lat() is None:
        continue
    points.append((school.lon(), school.lat()))


idx = rtree.index.Rtree()
for i, p in enumerate(points):
    idx.insert(i, p + p, p)

p = (way[-1].x, way[-1].y)
five_nearest = list(idx.nearest(p, 5, objects='raw'))


orig_point = Point(way[-1])
lines = [LineString([orig_point, Point(dest)]) for dest in points if dest in five_nearest]


geod = Geod(ellps="WGS84")
lines_length = [geod.geometry_length(line) for line in lines]
print(f'{lines_length} расстояния от точки до 5 ближайших школ в Бангладеше')

[2284559.9538974287, 2266832.2552604834, 2281156.469616435, 2281167.7889800607, 2267737.838541723] расстояния от точки до 5 ближайших школ в Бангладеше


sum_length = 0
for i in lines_length:
    sum_length += i
print(f'{sum_length} - сумма расстояний в метрах от точки до 5 ближайших школ в Бангладеше')

11381454.306296133 - сумма расстояний в метрах от точки до 5 ближайших школ в Бангладеше


import overpy
import shapely.geometry as geometry
from shapely.ops import linemerge, unary_union, polygonize


all_polygons = []
for area in areas:
    query = """[out:json][timeout:25];rel(""" + str(area[0]) +""");out geom;>;out skel qt; """ 
    api = overpy.Overpass()
    result = api.query(query)
    lss = [] #convert ways to linstrings

    for ii_w, ways in enumerate(result.ways):
        ls_coords = []
        for node in ways.nodes:
            ls_coords.append((node.lon, node.lat)) # create a list of node coordinates
        lss.append(geometry.LineString(ls_coords)) # create a LineString from coords
        
    merged = linemerge([*lss]) # merge LineStrings
    borders = unary_union(merged) # linestrings to a MultiLineString
    polygons = list(polygonize(borders))
    all_polygons.append([area[0], area[1], geometry.MultiPolygon(polygons)])


south_and_west_for_countries = []
for multypolygon in all_polygons:
    if multypolygon[1] == "France - Luxembourg" or multypolygon[1] == "France - Andorra":
        continue
    multy_poly = multypolygon[2]
    x = multy_poly[0].exterior.coords.xy[0][0]
    y = multy_poly[0].exterior.coords.xy[1][0]
    south = [x, y]
    north = [x, y]
    for poly in multy_poly:
        xs = poly.exterior.coords.xy[0]
        ys = poly.exterior.coords.xy[1]
        for x, y in zip(xs, ys):
            if y < south[1]:
                south = [x, y]
            if y > north[1]:
                north = [x, y]
    south_and_west_for_countries.append([south, north, multypolygon[0], multypolygon[1]])


railway_station = []
for country in south_and_west_for_countries:
    areaId = 3600000000 + country[2]
    overpass = Overpass()
    query = overpassQueryBuilder(area=areaId, elementType=['node'], selector='"railway"="station"', out='geom')
    result = overpass.query(query, timeout=900)
    count = 0
    if result.countElements() == 0:
        continue
    for node in result.nodes():
        geod = Geod(ellps="WGS84")
        south_line_length = geod.geometry_length(LineString([Point(country[0]), Point([node.lon(), node.lat()])]))
        north_line_length = geod.geometry_length(LineString([Point(country[1]), Point([node.lon(), node.lat()])]))
        if south_line_length < north_line_length:
            count += 1
    railway_station.append([count, country[2], country[3]])


bus_stop = []
for country in south_and_west_for_countries:
    areaId = 3600000000 + country[2]
    overpass = Overpass()
    query = overpassQueryBuilder(area=areaId, elementType=['node'], selector='"highway"="bus_stop"', out='geom')
    result = overpass.query(query, timeout=900)
    count = 0
    if result.countElements() == 0:
        continue
    for node in result.nodes():
        geod = Geod(ellps="WGS84")
        south_line_length = geod.geometry_length(LineString([Point(country[0]), Point([node.lon(), node.lat()])]))
        north_line_length = geod.geometry_length(LineString([Point(country[1]), Point([node.lon(), node.lat()])]))
        if south_line_length < north_line_length:
            count += 1
    bus_stop.append([count, country[2], country[3]])


max_ratio = [railway_station[0][0] / bus_stop[0][0]] 
for i in railway_station:
    for j in bus_stop:
        if i[2] == j[2]:
            if j[0] == 0:
                continue
            ratio = i[0] / j[0]
            if ratio > max_ratio[0]:
                max_ratio = [ratio, i[0], j[0], i[1], i[2]]
print(f'В {max_ratio[4]} (id: {max_ratio[3]}) самое большое отношение числа железнодорожных станций к автобусным остановкам в южной части - {max_ratio[0]}')
cameroon_id = max_ratio[3]

В Cameroon (id: 192830) самое большое отношение числа железнодорожных станций к автобусным остановкам в южной части - 0.6181818181818182


areaId = 3600000000 + cameroon_id
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"amenity"="cinema"')
result = overpass.query(query, timeout=900)

points = []
for cinema in result.elements():
    if cinema.lon() is None or cinema.lat() is None:
        continue
    points.append((cinema.lon(), cinema.lat()))


idx = rtree.index.Rtree()
for i, p in enumerate(points):
    idx.insert(i, p + p, p)

p = (way[-1].x, way[-1].y)
five_nearest = list(idx.nearest(p, 5, objects='raw'))

orig_point = Point(way[-1])
lines = [LineString([orig_point, Point(dest)]) for dest in points if dest in five_nearest]

geod = Geod(ellps="WGS84")
lines_length = [geod.geometry_length(line) for line in lines]


sum_length = 0
for i in lines_length:
    sum_length += i
print(f'{sum_length} - сумма расстояний в метрах от точки до 5 ближайших кинотеатров в Камеруне')

19612688.559784684 - сумма расстояний в метрах от точки до 5 ближайших кинотеатров в Камеруне


import json
import requests
import warnings

import folium
import pyproj
from shapely import geometry
from shapely.geometry import Point, LineString, Polygon, MultiPolygon

from openrouteservice import client


warnings.filterwarnings("ignore")


url = 'https://data.cityofnewyork.us/api/views/a3nt-yts4/rows.json?accessType=DOWNLOAD'

def CreateBufferPolygon(point_in, resolution=10, radius=10):    

    sr_wgs = pyproj.Proj(init='epsg:4326')
    sr_utm = pyproj.Proj(init='epsg:32632') # WGS84 UTM32N
    point_in_proj = pyproj.transform(sr_wgs, sr_utm, *point_in) # unpack list to arguments
    point_buffer_proj = Point(point_in_proj).buffer(radius, resolution=resolution) # 10 m buffer
    
    # Iterate over all points in buffer and build polygon
    poly_wgs = []
    for point in point_buffer_proj.exterior.coords:
        poly_wgs.append(pyproj.transform(sr_utm, sr_wgs, *point)) # Transform back to WGS84

    return poly_wgs


# Set up the fundamentals
api_key = '5b3ce3597851110001cf624857c34d71aefe42ed95ec0fb54e33a8ea' # Individual api key
clnt = client.Client(key=api_key) # Create client with api key
school_json = requests.get(url).json() # Get data as JSON

map_params = {'location':([40.703165, -73.987598]),
              'zoom_start': 5}
map1 = folium.Map(**map_params)

# Populate a school buffer polygon list
schools_poly = []
for school_data in school_json['data']:
    school_coords = school_data[8].split()
    school_coords = [float(school_coords[1][1:]), float(school_coords[2][:-1])]
    folium.features.Marker(list(reversed(school_coords)),
                           popup='School point<br>{0}'.format(school_coords)).add_to(map1)
    
    # Create buffer polygons around construction sites with 150 m radius and low resolution
    school_poly_coords = CreateBufferPolygon(school_coords,
                                           resolution=2, # low resolution to keep polygons lean
                                           radius=150)
    schools_poly.append(school_poly_coords)
    
    school_poly_coords = [(y,x) for x,y in school_poly_coords] # Reverse coords for folium/Leaflet
    folium.vector_layers.Polygon(locations=school_poly_coords,
                                  color='#ffd699',
                                  fill_color='#ffd699',
                                  fill_opacity=0.2,
                                  weight=3).add_to(map1)


# GeoJSON style function
def style_function(color):
    return lambda feature: dict(color=color,
                              weight=3,
                              opacity=0.5)

# Create new map to start from scratch
map_params = {'location':([40.744696, -73.953993]),
              'zoom_start': 12}
map2 = folium.Map(**map_params)

request_params = {'coordinates': [[-73.963793, 40.777111],
                                 [-73.989596, 40.703213]],
                'format_out': 'geojson',
                'profile': 'driving-car',
                'preference': 'shortest',
                'instructions': 'false',}
route_normal = clnt.directions(**request_params)
folium.features.GeoJson(data=route_normal,
                        name='Route with school',
                        style_function=style_function('#FF0000'),
                        overlay=True).add_to(map2)

school_buffer_poly = []
for school_poly in schools_poly:
    poly = Polygon(school_poly)
    school_buffer_poly.append(poly)


# Add the school polygons to the request parameters
request_params['options'] = {'avoid_polygons': geometry.mapping(MultiPolygon(school_buffer_poly))}
route_detour = clnt.directions(**request_params)

#Buffer route with 0.009 degrees (really, just too lazy to project again...)
route_buffer = LineString(route_detour['features'][0]['geometry']['coordinates']).buffer(0.009)
folium.features.GeoJson(data=geometry.mapping(route_buffer),
                        name='Route Buffer',
                        style_function=style_function('#FFFF00'),
                        overlay=True).add_to(map2)

for school_poly in schools_poly:
    poly = Polygon(school_poly)
    if route_buffer.intersects(poly):
        school_poly_coords = [(y,x) for x,y in school_poly] # Reverse coords for folium/Leaflet
        folium.vector_layers.Polygon(locations=school_poly_coords,
                                      color='#ffd699',
                                      fill_color='#ffd699',
                                      fill_opacity=0.2,
                                      weight=3).add_to(map2)
        folium.features.Marker(list(reversed(poly.centroid.coords[0]))).add_to(map2)

folium.features.GeoJson(data=route_detour,
                        name='Route without schools',
                        style_function=style_function('#00FF00'),
                        overlay=True).add_to(map2)

<folium.features.GeoJson at 0x7fe054a59dd0>


map2


# GeoJSON style function
def style_function(color):
    return lambda feature: dict(color=color,
                              weight=3,
                              opacity=0.5)

# Create new map to start from scratch
map_params = {'location':([40.744696, -73.953993]),
              'zoom_start': 12}
map3 = folium.Map(**map_params)

request_params = {'coordinates': [[-73.963793, 40.777111],
                                 [-73.989596, 40.703213]],
                'format_out': 'geojson',
                'profile': 'driving-car',
                'preference': 'shortest',
                'instructions': 'false',}
route_normal = clnt.directions(**request_params)
folium.features.GeoJson(data=route_normal,
                        name='Route without construction sites',
                        style_function=style_function('#FF0000'),
                        overlay=True).add_to(map3)

# Plot which schools fall into the buffer Polygon
school_buffer_poly = []
for school_poly in schools_poly:
    poly = Polygon(school_poly)
    folium.features.Marker(list(reversed(poly.centroid.coords[0]))).add_to(map3)
    school_poly_coords = [(y,x) for x,y in school_poly] # Reverse coords for folium/Leaflet
    folium.vector_layers.Polygon(locations=school_poly_coords,
                                  color='#ffd699',
                                  fill_color='#ffd699',
                                  fill_opacity=0.2,
                                  weight=3).add_to(map3)
    school_buffer_poly.append(poly)


# Add the school polygons to the request parameters
request_params['options'] = {'avoid_polygons': geometry.mapping(MultiPolygon(school_buffer_poly))}
route_detour = clnt.directions(**request_params)

folium.features.GeoJson(data=route_detour,
                        name='Route with construction sites',
                        style_function=style_function('#00FF00'),
                        overlay=True).add_to(map3)

# map3.add_child(folium.map.LayerControl())

<folium.features.GeoJson at 0x7fe04c482fd0>


route_buffer = LineString(route_detour['features'][0]['geometry']['coordinates']).buffer(0.009)
folium.features.GeoJson(data=geometry.mapping(route_buffer),
                        name='Route Buffer',
                        style_function=style_function('#FFFF00'),
                        overlay=True).add_to(map3)

<folium.features.GeoJson at 0x7fe0326cb090>


map3

Анализ пространственных данных. Домашнее задание №2¶

Задание №1. Горячая точка (алгоритм - 10 баллов, визуализация - 10 баллов).¶

Задание №2. Качество жизни (20 баллов).¶

Задание 2.1¶

Задание 2.2¶

Задание 2.3¶

Задание №3. Поездка по Нью-Йорку (маршрут - 20 баллов, визуализация - 10 баллов).¶